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ABSTRACT 

Leapfrog integration has been the method of choice in N-body simulations owing to 
its low computational cost for a symplectic integrator with second order accuracy. We 
introduce a new leapfrog integrator that allows for variable timesteps for each particle 
in large N-body simulations. Tests with single particles in fixed potentials show that 
it behaves as a symplectic integrator. We then examine the results of both standard 
leapfrog and our temporally adaptive leapfrog on full N-body integrations of clusters 
and large scale structure establishing accuracy criteria for both methods. The adaptive 
method shows significant speed-ups over single step integrations — but the integrator 
no longer appears to be symplectic or, in the case of large scale structure simulations, 
accurate. This loss of accuracy appears to be caused by the way that the timestep 
is chosen, not by the integrator itself. We present a related integration technique 
that does retain sufficient accuracy. Although it is not symplectic, it is apparently 
better than previous implementations and is our current integrator of choice for large 
astrophysical simulations. We also note that the standard leapfrog difference equations 
used in cosmological N-body integrations in comoving coordinates are not symplectic. 
We derive an implementation of leapfrog that is in comoving canonical coordinates to 
correct for this deficiency. 

Subject headings: Methods: numerical 
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1. Introduction 

Over the past decade, spatially adaptive meth- 
ods have been developed to calculate gravita- 
tional forces in N-body simulations. These in- 
clude tree codes (Appel, 1985; Barnes and Hut, 
1986), fast multipole methods (Greengard, 1987), 
and adaptive grids (Couchman, 1991). As the 
number of particles in an N-body simulation 
grows, so do the density contrasts. Hierarchi- 
cal methods can follow extremely large dynamic 
ranges in densities at modest additional cost per 
force evaluation. However, large ranges in den- 
sities also imply a large range in time scales 
(oc 1/yJ density) . If we take the final state of a 
simulation and weight the computational work 
done on particles not uniformly but inversely 
with their natural timesteps, we find a potential 
gain of ~ 50. Temporal adaptivity is one of the 
last algorithmic areas where we can target an or- 
der of magnitude improvement. Hence, we seek a 
hierarchical integrator, i.e., a method such that 
particles are on adjustable individual timesteps. 

The most commonly used time integration 
scheme for N-body simulations is the leapfrog 
method. Leapfrog has several advantages over 
other methods. 1) For second order accuracy 
only one force evaluation and one copy of the 
physical state of the system is required. This 
is particularly beneficial for N-body simulations 
where the cost of a force evaluation is very expen- 
sive. 2) The force field in an N-body simulation is 
not very smooth, so higher order does not neces- 
sarily mean higher accuracy. 3) It is a symplec- 
tic integrator, i.e., it preserves properties spe- 
cific to Hamiltonian systems. See Channell and 
Scovel (1990) for a review of symplectic integra- 
tors. Gravitational N-body systems are Hamilto- 
nian, and therefore they should benefit from the 
use of an integrator that conserves phase space 
volume and has no spurious dissipation. This 
could be especially important in self-gravitating 
systems where a dissipation time scale that is 
linked to the dynamical time can lead to a run- 



away to spuriously large densities. Leapfrog is a 
second order symplectic integrator requiring only 
one costly force evaluation per timestep and only 
one copy of the physical state of the system. 
These properties are so desirable that we con- 
centrate on making an adaptive leapfrog. 

Hierarchical leapfrogs have been used before 
(Porter 1985; Ewell 1988; Hernquist and Katz 
1989), but they were not symplectic (§4). There 
exist symplectic integrators with individual but 
fixed timesteps for either particles or modes (Saha 
& Tremaine 1994; Skeel & Biesiadecki 1994; 
Lee, Duncan & Levison 1997). Hut, Makino & 
McMillan (1995) proposed an iterative scheme 
for choosing timesteps for a single particle in a 
rapidly varying potential, but each iteration in- 
volves a force evaluation that is prohibitively ex- 
pensive for large N simulations (§2). 

Section 2 describes the theory behind a hierar- 
chical integrator. Section 3 describes tests of this 
integrator on a single particle in a potential, and 
Section 4 will present tests of the integrator on 
full N-body systems. We discuss the implications 
of these results in Section 5. 

2. Symplectic Integrators 

A symplectic integrator is an exact solution 
to a discrete Hamiltonian system that is close to 
the continuum Hamiltonian of interest. There- 
fore, it preserves all the Poincare invariants and 
places stringent conditions on the global geom- 
etry of the dynamics. An obvious example is 
total energy conservation in a system with a 
time-independent Hamiltonian or the conserva- 
tion of angular momentum in axisymmetric sys- 
tems (Zhang & Skeel 1995). A symplectic in- 
tegrator will exactly conserve the energy in the 
discrete Hamiltonian that is an approximation to 
the true energy of the system. This approximate 
energy oscillates about the true energy without 
any numerical dissipation. 

The difference between the discrete and con- 
tinuum Hamiltonians can be viewed clS Si small 
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Fig. 1. — A comparison of symplectic and non- 
symplectic integrators is made. The squares are a 
second order leapfrog integrator; the crosses are a 
4th order Runge-Kutta integrator with the same 
timestep, and the solid line is the exact solution. 

perturbation given by the truncation error of the 
integrator. The error is a Hamiltonian! If the 
error Hamiltonian is a sufficiently small pertur- 
bation, then the KAM theorem (Arnold 1978) 
guarantees that the invariant curves destroyed 
are a set of finite measure. In other words, al- 
most all orbits that are stable in the real sys- 
tem will continue to be stable in the numerical 
system. An illustration of these advantages is 
shown in Figure ||. Here the radial velocity, v r , 
is plotted against the radius, r, for an ellipticity, 
e = 0.5 Kepler orbit using a leapfrog integrator 
and using a fourth order Runge-Kutta integra- 
tor. In each integration, approximately 24 steps 
were taken per orbit, and the integrations ran for 
16 orbits. Note how the leapfrog integrator oscil- 
lates about the true solution but always remains 
on a one dimensional surface. This indicates that 
it is indeed conserving an energy-like quantity, 
i.e. having the orbit constrained to a one dimen- 



sional surface shows the existence of an isolat- 
ing integral of motion. On the other hand, the 
Runge-Kutta orbit slowly becomes more circular. 
The poor performance of the Runge-Kutta inte- 
grator is remarkable given that it is a fourth or- 
der integrator and uses four times as many force 
evaluations as the leapfrog integrator. Also note 
the large wiggles in the leapfrog integration at 
apoapse. These are indicative of the proximity 
of resonant islands that would lead to a instabil- 
ity for larger timesteps. 

The leapfrog integrator can be written as 

r n+i/2 = r„ + ^rv n , 
Vn+l = v„ + Ta(r n+1 / 2 ) , 
Tn+l = r n+1 / 2 + 5 TV n+l) 

where r is the position vector of a particle, v is 
the velocity, a is the acceleration, and r is the 
timestep. When several of these steps are put 
together, the two position updates can be com- 
bined to a single update r n _]/ 2 to r n+1 / 2) and the 
resulting alternation between updating r and v 
gives leapfrog its name. The symplectic nature of 
leapfrog can be seen by noting that the position 
update is equivalent to evolving the system ex- 
actly under the Hamiltonian Hp = \v 2 , and the 
velocity update is equivalent to evolving the sys- 
tem under the Hamiltonian Hr = V(r), where 
V(r) is the potential generating the accelerations. 
We will call the operator that evolves the system 
under Hp the "drift" operator, D, and the opera- 
tor that evolves the system under Hk the "kick" 
operator, It can be shown that (see Saha 

and Tremaine, 1992 for details) the combination 
of operators D(r /2)K{t)D{t /2) evolves the sys- 
tem under a Hamiltonian 

Hn = Hd + Hk + H err = -v 2 + V(r) + H err , 

where H err is of order r 2 . The existence of this 
surrogate Hamiltonian ensures that the leapfrog 

2 If one expresses Hamilton's equations as z = {z,H}, 
where z are the phase space coordinates and { , } are 
Poisson brackets, then formally D(t)z = exp(r{z, Hn}), 
and K(r)z = exp(r{z, Hk}). 
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is symplectic and second order. Since the Hamil- 
tonian is symmetric with respect to Hp and Hk, 
the combination of operators K(t/2)D(t)K(t/2) 
is also a symplectic second order integrator. Ex- 
plicitly, this is 

v n+ i/2 = v n + ^ra(r„), 
r n +i =r n + rv n+1 / 2 , 
v n +i = v n+1/2 + ira(r n+ i). 

Higher order symplectic integrators can be con- 
structed from combinations of leapfrog steps (Yoshida 
1990). Each Nth. order leapfrog integrator re- 
quires N — 1 force evaluations and only one copy 
of the physical state of the system. 

Unfortunately, constructing a variable step- 
size method by choosing a new timestep after 
each leapfrog step gives very disappointing re- 
sults (Calvo and Sanz-Serna, 1993). Some sim- 
ple schemes can be shown to have no rigorous 
stability criterion for any step size (Skeel 1993). 
Several explanations have been given for this be- 
havior. The simplest is to note that a variable 
step integrator is evolving a dynamic system with 
state variables (r, v, t). The projection onto 
the phase space coordinates (r, v) cannot be de- 
scribed by a Hamiltonian. Skeel and Gear (1992) 
consider a one step symplectic operator of the 
form 

(r„+i,v n+ i) = F(r n ,v n ,r(r n ,v n )), 

and show that if F is symplectic for a constant 
step size, it will not in general be symplectic for 
variable step size. Another way to see this prob- 
lem is to notice that the time reversibility has 
been broken: if we step forward in time and then 
step backward, we do not end up at the same 
point because of the change in timestep. Since 
symplectic implies time reversible, such an oper- 
ator is not symplectic. Integrators that are time 
reversible are referred to as reflexive by Kahan 
(1993; see also Kahan and Li 1997) who argues 
that reflexivity is the key to the robust properties 
of "updating formulae" rather than symplectic- 
ity. 



A strategy to make leapfrog reflexive are clear 
from a cursory look at its evolution operators. 
The operation D(t/2)K(t)D(t/2) is reversible 
if we use an "select" operator to choose the 
timestep, S, such that the time reversibility is 
retained. Hut, Makino, and McMillan (1995) 
achieve this by using an implicit definition of the 
timestep: 

T = 7}i T ( r n,v n ) +r(r n+ i,v n+ i)]. 

So, the beginning and end of the step are re- 
quired to agree on the timestep. One can solve 
for r iteratively, and very good results are ob- 
tained even with only one or two iterations. (This 
result is also clear from Kahan 1993). This itera- 
tive scheme poses several problems when applied 
to a large N-body simulation. It requires backing 
up timesteps, throwing away expensive force cal- 
culations, using auxiliary storage and must spec- 
ify a method for synchronizing the particles for 
mutual force evaluations. 

However, we can also see an alternative means 
to restore reflexivity. Let us choose a select op- 
erator S that commutes with K, so that DSKD 
is equivalent to DKSD. Since K only changes 
the velocities and not the positions, an S op- 
erator that depends entirely on positions satis- 
fies the commutation requirement. As we shall 
show below, this is not strictly time reversible. 
However, since the operators read the same for- 
ward and backward, we will refer to DSKD 
as a "palindromic" integrator. Synchronization 
can be maintained by only choosing timesteps 
that are a power-of-two subdivision of the largest 
timestep, r s . That is, 

- Il 

where Tj is the timestep of a given particle, and 
rii is an integer (See Hernquist and Katz, 1989). 
Combining this synchronization procedure with 
the DSKD evolution operation we have the fol- 
lowing (recursive) algorithm for a timestep. 1) 
Drift the particles forward Tj/2. 2) Apply the 
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select operator. If it accepts this step, then fin- 
ish the step with K{ri) and D(n/2). Other- 
wise, drift the particles back n/2 and take two 
timesteps using Tj/2. The algorithm can be ex- 
pressed recursively in the following pseudo-code: 
Timestep(r) 

{ 

Drift (r) ; 
Select (r) ; 

if T ~ Train 
{ 

Kick(r) ; 
Drift (r) ; 

} 

else 

{ 

Drift (-r) ; 
Timestep(r/2) ; 
Kick(r) ; 
Timestep(r/2) ; 

} 

} 

Here, Drift (r) applies D (r) to all particles, 
Kick(r) applies K{t) to particles with timestep 
r and Select () decides to which timestep the 
particles belong, with the side effect of finding 
the minimum timestep, r m j n . We will refer to 
the traditional method of choosing a timestep at 
the beginning of the integration step as SDKD. 

As an example of how this procedure works, 
consider Figure ||. In this diagram we con- 
sider the case where there are particles on three 
separate timesteps. The arcs represent drift- 
ing the particles, and the vertical slashes repre- 
sent either an select or kick of particles at that 
timestep. Starting at point 0, all particles are 
drifted to point 1 and evaluated for whether they 
are on the largest timestep. Then all particles 
are drifted to point 2 and those particles not on 
the largest timestep are evaluated for whether 
they are on the middle timestep. Then all parti- 




Fig. 2. — Diagram of particles on three separate 
timesteps, where arcs represent "drifts", and ver- 
tical lines represent "kicks". For DSKD, the or- 
der of flow throw the diagram is 0, 1, 2, 3, 2, 4, 
1, 5, 6, 5, 7, 8. 

cles are drifted to point 3 where those particles 
on the smallest timestep are kicked. All parti- 
cles are then drifted to point 2 where particles 
on the middle timestep are kicked, then to point 
4 where particles on the smallest timestep are 
kicked, then to point 1 where particles on the 
largest timestep are kicked. Then all particles 
are drifted to point 5 where again all particles not 
on the largest timestep are evaluated for whether 
they are on the middle timestep. Then the ap- 
propriate particles are kicked at points 6, 5, and 
7. Finally all particles are drifted to point 8 and 
a single large timestep is complete. 

There are several things to note about the 
above algorithm. One is that it is not unique. It 
uses a "top down" approach by trying the largest 
timestep and reducing it until it is deemed ap- 
propriate. A "bottom up" scheme can be imple- 
mented where a smallest timestep is tried first. 
Secondly, this algorithm is not exactly time re- 
versible. Scenarios can easily be constructed 
where a forward step followed by a backward step 
will not come back to the initial conditions. An 
example is shown in Figure [| The circle de- 
lineates the boundary between region 1, where 
one timestep is needed, and region 2, where a 
timestep of 1/2 the timestep in region 1 is needed. 
If we start a particle in region 2, it is drifted 
forward, found to still be in region 2, so two 
timesteps are taken. If we start the particle in 
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Fig. 3. — An example of a non-time-symmetric 
scenario is shown. If the step starts in region 2, 
two steps are taken. If the step starts in region 
1, only one step is taken. 



DSKD, 77= 0,1 
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region 1 and go backward, it is possible that the 
drift leaves the particle still in region 1, so only 
one timestep is taken going backward, and we do 
not arrive at the point from which we started. 
Finally, since several selects may be done for ev- 
ery kick, the efficiency of this algorithm depends 
upon Select being able to quickly decide an ap- 
propriate timestep given the particle positions. A 
suitable criterion that depends only on the posi- 
tions is one based upon the local dynamical time 
td = 1/ 1 \JGp. With this criterion, Select will 
pick the largest timestep r such that 



where rj is a constant to be determined based on 
stability and accuracy requirements. With this 
criterion, region 2 in Figure | is a high density 
region while region 1 has low density. Similar 
problems to this loss of reflexivity can occur in 
the iterative schemes (Hut, Makino, & McMillan 
1995; Kahan 1993) if the starting guess of an 
iteration is preconditioned by the history of the 
particle. That is, the loss of reflexivity depends 
as much on the choice of a "top down" approach 
as on the Select operator. 



Fig. 4. — Relative change in total energy is plot- 
ted as a function of time in units of the orbital 
period for several integrators. Curves are plot- 
ted for DSKD with 77 = 0.1 and rj = 0.03, and 
SDKD with r, = 0.03. 

3. Single Particle Tests 

Our first test is the motion of a test particle 
in the one dimensional effective potential of the 
Kepler problem. If the integrator performs well, 
these tests will guide our choice of 77 for later 
applications. The timestep criterion in equa- 
tion |l] is not straight-forward to implement in 
this case, since the density is zero anywhere out- 
side the central source. Therefore, in this test 
we determine the timestep by the enclosed den- 
sity, p e = 3M/47rr 3 where M is the mass of the 
central object. 

Figure || shows the relative change in total en- 
ergy as function of time for several integrators. 
The orbit has eccentricity e = 0.5, and is evolved 
for about 100 periods. Results for two values of 
rj for the DSKD method are plotted as well as a 
result for SDKD. A fixed step integration was 
also run, and its energy nearly coincides with the 
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rj = .03 DSKD method. Prom the energy plot, it 
appears that the DSKD method with r] = 0.03 
is indeed symplectic. In contrast, the SDKD 
method with the same timestep criterion experi- 
ences a secular drift in the energy. The savings 
in force evaluations is also significant. A fixed 
timestep integration with similar energy errors 
required 50000 force evaluations to accomplish 
100 orbits while the adaptive timestep integra- 
tion required less than 16000, a savings of over 
a factor of 3 for this moderate eccentricity orbit. 
Also, it can be seen from the energy plot that 
i] ps 0.03 is needed for the adaptive algorithm to 
be stable. This is because a stable orbit, having 
one integral of motion for each degree of free- 
dom, should be quasiperiodic, which the r\ = 0.03 
line appears to be.^| For r\ larger than this, the 
method does not appear to be symplectic, but 
the changes in energy appear to be due to insta- 
bility, (i.e. integrals being destroyed) rather than 
a secular buildup of truncation error. 

A more realistic test of the integrator is the 
motion of a single particle in the potential of an 
isothermal sphere. As well as having a density 
defined everywhere, allowing the use of the lo- 
cal density in the timestep criterion, the poten- 
tial is similar to that of globular clusters, galax- 
ies and clusters of galaxies in our large N-body 
simulations. Figure || shows the energy conser- 
vation for the SDKD and DSKD methods for 
a particle in the potential of a singular isother- 
mal sphere. The ratio of apocenter to pericen- 
ter is 3.2:1. The density used for choosing the 
timestep is the local density of the isothermal 
sphere. Again the SDKD shows a secular drift 
in the energy, while the DSKD has the charac- 
teristics of a symplectic method. Also note that 
in the case of an isothermal sphere, the method 
is stable for r] = 0.1. 

3 Remember that if there is one integral for each degree of 
freedom, the Hamiltonian can be expressed as H = H(J), 
and any function of the phase space coordinates, f(J,9) 
can be expressed as f(t) — /(J,ut) where J and u> = 
dH/dJ are constants. 
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Fig. 5. — Relative change in total energy is plot- 
ted as a function of time in units of the orbital 
period for several integrators in an isothermal 
potential. Curves are plotted for DSKD with 
rj = 0.1 and 77 = 0.03, SDKD with rj = 0.03. 

4. N-Body Tests 

For our purposes, the real proof of a given in- 
tegration method is how well it performs in an 
actual N-body code. In this section we test the 
above methods in two complementary systems. 
The first is a King model — a model which, in 
the absence of collisional relaxation, should re- 
main static. The second is a cosmological sim- 
ulation of a universe dominated by cold dark 
matter (CDM). This case is very dynamic, with 
structure evolving on all scales. All the simu- 
lations below were performed with PKDGRAV 
(Dikaiakos & Stadel, 1996; Stadel & Quinn, in 
preparation), a parallel code that uses a binary 
tree with a Barnes and Hut (1986) style opening 
criterion to make the gravity calculation order 
N\og{N). To mitigate the effect of force errors, 
we use an opening criterion 9 = 0.55, and accel- 
erations from cells are expanded to hexadecapole 
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Fig. 6. — The spherically averaged density in 
units of the central density is plotted against ra- 
dius in units of the core radius for the end state of 
several simulations. The solid line is the density 
profile of the initial King model. All integrations 
were done with a fixed step leapfrog integrator 
with the timestep in units of the central dynam- 
ical time given in the legend. 

order in all the following simulations. 

Before we test the multistep model in a full 
N-body simulation, let us explore the more basic 
question of what single fixed timestep is appro- 
priate for an astronomically interesting N-body 
model. For our tests we will use a W = 9 King 
model realized with 100,000 particles evolved for 
100 central dynamical times. This is representa- 
tive of a galaxy halo over its lifetime in a typical 
cosmological N-body simulation. 

Energy conservation is a typical measure of 
the quality of a simulation, but more appropri- 
ate measures should involve the convergence of 
scientifically interesting quantities. For the King 
model tests we will use the radial density profile 
which should remain constant in a collisionless 
simulation. In Figure ^ we compare the density 
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Fig. 7. — The relative change in the total energy 
as a function of time for the same simulations as 
in Figure |(| The time is in units of the central 
dynamical time. 

profile of the initial King model with that at the 
end of a simulation with various timestep sizes 
evolved over 100 central crossing times. All mod- 
els were integrated with fixed step leapfrog using 
timesteps of either 1.0, 0.3, 0.1, or 0.03 times the 
central dynamical time. From the figure, we see 
that t < O.ltd is needed to maintain the cen- 
tral density of the King model. Figure |7] shows 
the evolution of the total energy for these same 
simulations. Note that an energy conservation of 
better than 3% is needed to preserve the density 
profile of this King model. 

In order to use a density criterion to deter- 
mine timesteps in a general N-body simulation, 
it is necessary to have a method for calculating 
the local density for every particle in an arbitrary 
distribution. We use a Smooth Particle Hydro- 
dynamics estimate of the density by smoothing 
over a fixed number of nearest neighbors with a 
cubic spline kernel. The smoothing length, h, 
is adaptive and is set so that there are exactly 
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Fig. 8. — The spherically averaged density is 
plotted against radius as in Figure |(| The in- 
tegrations were done using a multistep leapfrog 
with a density timestep criterion. All integra- 
tions except the one indicated used a DSKD 
style timestepping. 

64 particles within 2h. The kernel is made sym- 
metric using the "gather-scatter" algorithm as 
described in Hernquist and Katz (1989). 

We show the radial density profiles for adap- 
tive timestep integrations in Figure |8L From 
the figure, we see that, as in the single particle 
case, T) < 0.03 is needed to integrate this model 
with reasonable accuracy. The DSKD integra- 
tion with rj = 0.03 needed a factor of 6 less force 
evaluations than the fixed step integration with 
r/td = 0.03, or a factor of 2 less than the fixed 
step, r/td = 0.1 integration. As shown in Figure 
^, the energy conservation of the adaptive scheme 
is comparable to the fixed step integrator. Also 
shown in the figures are results for the SDKD 
integrator. It appears that at the level of accu- 
racy needed to integrate this King model for 100 
dynamical times, the secular drifts introduced by 
the SDKD integrator do not significantly affect 



Fig. 9. — The relative change in the total energy 
as a function of time for the same simulations as 
in Figure [8|. The time is in units of the central 
dynamical time. 

the integration. 

Before testing timestepping in a cosmologi- 
cal simulation we first note that the standard 
leapfrog difference equations used in cosmolog- 
ical simulations with comoving coordinates are 
not done in canonical coordinates and are appar- 
ently not symplectic. However, as shown in the 
appendix, a symplectic integrator can easily be 
derived, for which the "drift" and "kick" opera- 
tors are 

r t+T dt 

D(r)^v' t+T =v' t+V ' / - 2 

Jt a z 

r t+T dt 

k (r) = p' t+ . = p' t - vy / -, 

Jt a 

where the variables are as defined in the ap- 
pendix. We use this integrator in the tests that 
follow. This integrator has another advantage 
over the standard leapfrog difference equations in 
that it can be used to implement both the DKD 
and KDK form of leapfrog. The standard equa- 
tions can only be derived for DKD. The KDK 
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form of leapfrog can have advantages over DKD 
(see the discussion of costs in §5). 

Energy conservation may be a an even poorer 
criterion for determining the accuracy of a cos- 
mological integration. This owes to the possi- 
bility of making small errors in calculating the 
energy of the uniform background that can dom- 
inate the total error budget and mask any prob- 
lems in the structure of interest. Nonetheless, it 
pays to examine the evolution of the energy. It 
was transients in the energy at the beginning of 
the cosmological simulations that made us sus- 
pect that the integration was not in canonical 
coordinates, leading to the results in Appendix 
A. But, these large transients in energy had little 
effect on the three dimensional structure, under- 
scoring our point that it may be a poor diagnos- 
tic. Small changes in the way that two simula- 
tions handle global energy may be an interesting 
clue to other aspects of the integrator. 

Since there is no analytic model for the full 
three dimensional formation of non-linear cos- 
mological structure, we have to use another sim- 
ulation to determine the ground truth. There- 
fore, we test the goodness of a given timestep- 
ping algorithm by examining its convergence to a 
simulation with a large number of timesteps for 
all particles. This convergence will be assessed 
based on properties of the groups that are formed 
and the internal structure of the largest group. 

The fiducial simulation is of a standard CDM 
dominated universe in a periodic box of size 
22.2 Mpc (h = 0.5) on a side realized on a grid 
of 32 3 particles. This implies a particle mass of 
2.3 x 10 10 Mq. A spline softening is used with 
a softening length of 20 kpc. The initial condi- 
tions were imposed on the particles at a redshift 
of 49 using the Zel'dovich approximation. The 
small size of this box means it does not accu- 
rately represent large scale structure, but it does 
guarantee significant non-linear evolution and a 
stringent test for the time integrator. 

The first convergence test is the density profile 
at the end of the simulation of the largest halo as 




Fig. 10. — The density averaged in spherical 
shells is plotted as a function of radius for the 
largest halo in the cosmological simulation. 

identified using a friends-of-friends grouping al- 
gorithm. The object in question has a total mass 
of 1.23 x 10 14 M Q within the virial radius and a 
maximum circular velocity of 639 kms -1 . This 
test looks at how the timestep algorithm affects 
the structure of well resolved objects in the sim- 
ulation. Since structure forms hierarchically, the 
objects in the largest clusters will have experi- 
enced a complicated history with several dynam- 
ical times in previous generations of structure. 
Any effect that accumulates with dynamical time 
is most likely to betray itself in the center of the 
richest cluster. 

Comparisons of the density profile of this ob- 
ject in various integrations are shown in Figure 
10. The fixed step run used 4000 steps in a Hub- 
ble time. The DSKD and SDKD runs used 
rj = 0.03. As might be expected, the timestep 
criterion used for the King model simulation per- 
forms well in accurately modeling the density 
profile of this cluster-size object. Looser crite- 
ria (rj > 0.1) fail to capture the central density 
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Fig. 11. — me cumulative mass traction m 
groups of a given particle number or higher is 
plotted. 

of the object, just as in the static King model. 

Our second criterion is the cumulative mass 
fraction in groups, also as identified by the friends- 
of- friends algorithm with a linking length of 0.26 
times the mean interparticle separation. This 
linking length corresponds to an enclosed density 
equal to the virial density assuming an isother- 
mal sphere. Here the motivation is to see the 
effect of the timestep algorithm on the small- 
est objects we can resolve. The cumulative mass 



fraction of halos is plotted in Figure 11. 

For the smallest objects, the density criterion 
does not do well. There is a significant decrease 
in the number of objects formed with less than 
one hundred particles in the simulation using the 
density timestep criterion. It may not be surpris- 
ing that particles residing in halos with compa- 
rable or fewer particles than that used by the 
smoothing kernel are unable to accurately esti- 
mate their local density. On the other hand, a 
simulation using a (non-time symmetric) multi- 
stepping algorithm based on the acceleration of 



the particles reproduces the numbers of objects 
seen in the fixed step simulation all the way down 
to 8 particles. (The discrepancy seen at about 
500 particles per group is due to a merger occur- 
ring at the end of the simulation.) 

5. Discussion 

First, we examine our results in terms of the 
number of timesteps needed to accomplish a 
given simulation with reasonable accuracy. We 
will determine the number of fixed steps before 
looking at the gain from hierarchical timestep- 
ping. For an equilibrium model where we wish 
to maintain the overall density profile, the an- 
swer is straightforward: the timestep should be 
set to0m/^/Gp 

max 5 where p m ax is the maximum 
density. For a cosmological simulation, more as- 
sumptions need to be made. Consider a simula- 
tion in which the largest halo has a circular veloc- 
ity, v c , that is roughly constant with radius down 
to some core radius, rQ. If we wish to accurately 
model this halo, then our timestep criterion must 
be 

r] Tyro / 47T 

vr 



At 



VGpo v c 

where po is the (roughly constant) central den- 
sity, and we have assumed spherical symmetry. 
Now consider a simulation with a halo with v c = 
1000 kms -1 where our goal is to resolve its struc- 
ture within 10 kpc. In units of the Hubble time, 
1/Hq, we have 



H At ^ 1 x 10' 



V 



ro 



Vr. 



-1 



,10 kpc7 V1000 kms" 

For 77 = 0.03, this gives 33,000 steps per Hubble 
time. Previous estimates of the timestep con- 
clude that 6000(10 kpc/e) timesteps are needed 
per Hubble time, where e is the gravitational soft- 
ening length (Lake et al. 1995). This is consistent 
with our estimate if we assume that the gravita- 
tional softening is 1/5 the core radius we wish to 
resolve. 

To evaluate the performance of a multistep in- 
tegrator, we must determine the "fixed costs" re- 
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quired for any level of the timestep hierarchy as 
well as the costs in determining the timestep level 
assigned to each particle. The total cost of force 
evaluations is then added and amortized over the 
longest timestep. For example, a tree code re- 
quires that the tree is built for all the particles 
regardless of the number of particles requiring 
force evaluations. Indeed, for a standard particle 
mesh code, the forces at every grid point are ei- 
ther calculated or not. The only part of the force 
evaluation that is not part of the fixed cost would 
be the trivial interpolation of the force onto par- 
ticles. So, essentially all the costs are fixed and 
no gain is possible. 

First, we consider the costs in the SDKD 
scheme with our tree code, PKDGRAV, and as- 
sume that the cost of actually moving the parti- 
cles and evaluating the timestep is negligible. A 
single base step with r different timesteps sepa- 
rated by a factor of two, will have a cost of 

r 

(2-_l)C t + C / ^2 i - 1 iV: j , 
1=1 

where Ct is the cost to build the tree, Cf is 
the cost of calculating the force on a single par- 
ticle, and Ni is the number of particles on a 
timestep i of size 2 1_ * times the base step. For 
a fixed timestep simulation with all particles at 
the smallest timestep, the cost would be 

2 r -\C t + NC f ). 

Comparing these two costs, it is easy to see that if 
Ct dominates the computation, then in the limit 
of large r, the SDKD scheme would a factor of 
two more expensive than the single step calcu- 
lation. On the other hand, if Ct is negligible, 
then the maximum speedup would be a factor of 
2 r_1 if almost all the particles were on the largest 
timestep. For a given ratio, /, between Ct and 
the cost of evaluating all the forces, NCf, there is 
a maximum speedup of (/+l)/2/ for large r, and 
almost all of the particles on the largest timestep. 
Note that for an SKDK scheme this maximum 



speedup would be a factor of two larger since 
gravity evaluations for larger timesteps are syn- 
chronized with those of smaller timesteps, and 
the same tree can be used for both. For PKD- 
GRAV on the King model with 100,000 particles, 
the ratio / is about 0.024, giving a maximum 
speedup of about 22. 

For the symplectic scheme with a density cri- 
terion, the cost of evaluating the timestep is non- 
negligible, and we must account for it. It also has 
a fixed cost part for building the tree, C s t, and 
a per particle cost, C s f. In terms of these costs 
the total cost of a base step is now 

r r 

(2 r -l){C t +C st )+C f ]T 2 i " 1 iV l +C s/ ^(tf-lJJVi. 

i=l i=l 

If the tree build costs and per particle costs for 
selecting the timestep and evaluating the forces 
are similar, then the symplectic scheme would 
be a factor of two more expensive than the non- 
symplectic scheme. For PKD GR A V running on 
a 100,000 particle King model, the per parti- 
cle cost of the density criterion is a factor of 19 
smaller than the per particle cost of a force eval- 
uation, and the tree build for a density tree is 
50% faster than for a gravity tree. The maximum 
speedup assuming an optimum particle distribu- 
tion is then about 13. 

For realistic particle distributions, the speed 
up can be much smaller, especially since a few 
particles (0.1% in the 100,000 particle King model) 
end up on timesteps smaller than the single step- 
ping timestep. The speedup of the multistep run 
over the single step run as calculated from the 
above formulae is about a factor of 4. Note that 
we have neglected some costs such as particle 
pushing and the inefficiency of calculating forces 
for a few particles on modern pipelined proces- 
sors. The actual factor in wall clock time is 2.7. 

For non-equilibrium situations such as cosmo- 
logical simulations, a multistep scheme offers ad- 
ditional speedups since it can adapt to the chang- 
ing time scales in the simulation. An exam- 
ple of the timestep evolution is shown in Figure 
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Fig. 12. — The cumulative fraction of particles on 
a given timestep or larger is plotted as a fraction 
of the age of the Universe in a cosmological simu- 
lation. The lowest line is for the longest timestep 
taken, and each upper line is for a factor of two 
smaller timestep. 

[T^ . Here, the fraction of particles at a partic- 
ular timestep or larger is plotted as a function 
of time. Note how almost all particles start out 
on a very small timestep at the beginning of the 
simulation when the Universe is very dense, and 
generally migrate to larger timesteps as the sim- 
ulation progresses. One could make adjustments 
for this trend in a single step simulation by trans- 
forming to an appropriate variable, e.g. expan- 
sion factor instead of time (as in Efstathiou et 
al. 1985), but as can be seen from the figure, a 
simple monotonic transformation would not cap- 
ture all the complexities of the changing time 
scales. Furthermore, one would need different 
transformations for different cosmological mod- 
els. Finally, to keep the integrations in canonical 
coordinates would require that the gravitational 
softening length is fixed in coordinates that were 
different from either comoving or physical (see 



Appendix A). 

From our pragmatic standpoint as users of N- 
body simulations to model complex phenomena, 
the gain from symplectic integrators is their abil- 
ity to tolerate greater truncation errors in numer- 
ical simulations of Hamiltonian systems. This, 
in turn, allows for longer timesteps and short- 
ens the computer time needed to complete the 
simulation. Similarly, individual timesteps in a 
particle simulation allows one to concentrate the 
computational effort on those particles with the 
shortest dynamical times that require it. Either 
of these is beneficial only if their advantages out- 
weigh their cost, whether in computational effort 
or algorithmic complexity. 

For a system like the Solar System that is sim- 
ulated for billions of dynamical times, a sym- 
plectic integrator has obvious advantages. Any 
dissipation introduced by truncation error has 
dire consequences that are easily observable — 
a planet may spiral into the Sun. The al- 
ternatives to symplectic integrators are either 
very high order integrators or extremely short 
timesteps, or both. Even a computationally ex- 
pensive symplectic integrator can be an overall 
winner against these alternatives. (Wisdom and 
Holman, 1991; Saha and Tremaine, 1992.) 

For other systems, such as galaxies or large 
scale structure, the situation is not so clear. 
These systems are simulated for at most a few 
hundred dynamical times, and drifts in conserved 
quantities due to truncation error will not be as 
noticeable. In these marginal cases, a symplectic 
integrator will only be beneficial if it is simple to 
implement and computationally cheap. In this 
paper we focus on a criterion based on the local 
density for this very reason: the algorithm for es- 
timating the local density is very fast compared 
to calculating gravity. Unfortunately, this cri- 
terion is not ideal, particularly for cosmological 
simulations. 

Is there a better criterion? If we use the kick- 
drift-kick form of leapfrog, one could imagine a 
timestep criterion that depends only on the ve- 
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locity, and therefore commutes with the drift op- 
erator. An example would be r < rje/v, where 
e is the force softening. Three objections can 
immediately be raised against such a criterion. 
First, it is not Galilean invariant. A particle mov- 
ing with respect to the coordinate system of the 
simulation will be on a different timestep than 
one that is at rest irrespective of the forces act- 
ing on the particle. Second, if the acceleration of 
a particle is nearly opposite to its velocity at the 
beginning of a timestep, then the timestep can 
be chosen such that the velocity is nearly zero 
at the middle of the timestep where the Select 
decision is made. Such particles will thereby be 
always put on a very long timestep. Third, the 
choice of e to set the scale assumes that the dy- 
namical times are set by two-body encounters. 
For a particle orbiting in a cluster, it is the clus- 
ter density rather than two-body encounters that 
set the dynamical time. On the other hand, in a 
cluster with significant substructure the timestep 
should be set by the encounters with this sub- 
structure. In this case, e/v might be an accurate 
estimate of the timestep. 

An obvious criterion is the acceleration, for 
example r < rj^e/a. This commutes with the 
Kick operator because Kick changes only the ve- 
locities, and the acceleration is only a function 
of the positions. So one could easily use it as 
the select operator in DSKD just like we use 
the density criterion. However, the acceleration 
is just the thing that we are trying to minimize 
calculating. Using it as a timestep criterion will 
cancel the computational savings we are trying 
to achieve with multistepping. For an SDKD 
method, this is not a problem since one could 
use the last calculation of the accelerations that 
were used to advance the velocities. This luxury 
is not available for the DSKD method, and one 
would have to evaluate the acceleration of a par- 
ticle many times during the Select phase. One 
could also imagine criteria based on higher order 
terms, such as the magnitude of the divergence 
of the acceleration, or the time derivative of ac- 



celeration. However, these tend to be even more 
expensive to calculate than the acceleration. 

Although our DSKD method appears to be 
symplectic when tested with single particle inte- 
grations, it is not the best integrator for general 
N-body problems. We simply do not have a suit- 
able Select operator that commutes with either 
the Drift or Kick operator, that is computation- 
ally cheap to evaluate compared to an accelera- 
tion, and that returns an accurate timestep in 
all cases. The local density criterion does well 
on the first count but only occasionally satisfies 
the second. This problem was foreshadowed by 
the inability to define the local density in the 
Kepler case and having to resort to a global cri- 
terion that used the enclosed density. When cal- 
culated in a cosmological simulation by averaging 
over the nearest N particles, it erased substruc- 
ture with < N particles. We have not exhausted 
the search for an ideal timestep criterion, and 
there yet may be something that fits these strin- 
gent requirements. Nevertheless, the criterion we 
have explored in this paper may be very useful 
in simulations where structure is defined by large 
numbers of particles, such as in the modeling of 
galactic structure. 

For more general cosmological N-body simula- 
tions we advocate either the SDKD or SKDK 
method using r)^/e/a and/or (if clusters have 
significant substructure) i](e/v) to choose the 
timestep. Using tests similar to those described 
for our symplectic integrator we determined that 
an appropriate choice for 7] is 0.3 (for Plummer 
softening 0.4) to insure stability of the integra- 
tion. Although this algorithm is not symplectic, 
it appears to give accurate results for quantities 
of interest in a cosmological N-body simulation. 

This work was supported in part by NASA 
HPCC/ESS grant NAG 5-2213 and by NASA As- 
trophysics Theory grant NAGW-2523. We also 
wish to acknowledge useful discussions with Scott 
Tremaine. 
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A. Symplectic Integrators in Comoving 
Coordinates 

The equations of motion in a comoving coor- 
dinate frame are traditionally presented as 



v' + 2H(t)V 



V 



r = v 



W = 47rG(p'-p / 6 ), 

where the primed quantities refer to the comov- 
ing coordinate frame, a is the expansion factor, 
H is Hubble's constant, and p' b is the mean back- 
ground density. These can be integrated using 
the difference equations 

r'n+l/2 = r' n + ^rv' n 

, 1 - H(t)T V>'(r' n+ 1/2)T 



V n+1 
r 



" v n l + H(t)r + a 3{l + H{t)T) 

'n+1 = r'n+l/2 + 2 rv 'n+l- 



(See Hockney and Eastwood, 1981.) However, 
these equations do not appear to be symplectic, 
i.e., there is no known generating function that 
produces this transformation. 

However, by making a suitable canonical trans- 
formation, one can easily derive an integrator 
that is symplectic. The Lagrangian for the par- 
ticle motion in the comoving frame is 

C = I( av ' + dr') 2 -<f>. 

With the (time dependent) generating function 
tp = (l/2)adr /2 , this transforms to 



C = ±a 2 v' 2 



if/fa, 



where 



b' = a<j> + ±aa 2 r' 2 : . 



(See Peebles, 1980.) 

Switching to the Hamiltonian formalism, the 
momentum canonical to r' is p' = a 2 v', and the 
Hamiltonian is 



2a 2 a 



Although this Hamiltonian is time-dependent^, it 
is separable, so the "drift" and "kick" operators 
are easily derived as: 

D(r) = r' t+T =r' t+P ' f +T % 

Jt a, 

r t+T dt 

K(r)E P ' i+r =p' ( -V'W -, 

Jt a 

where r is the step size. Note that it is assumed 
here that there is no explicit time dependence in 
(f)'; e.g. any softening of the potential must be 
constant in comoving coordinates. For standard 
cosmologies the integrals in the above operators 
can be easily evaluated. For example in a critical 
density universe, 

r*+ T dt 2 



Hn 



-1/2 



-1/2 



(t + r) 



and 



,1/2 



(t + r) 



,1/2 



(t) 



<- t+T dt _ 2 
it a H 

For non-flat matter dominated universes, it is 
convenient to use the parametric solutions (Pee- 
bles, 1980), 

a = A(l — cos 77), t = B{rj — sin 77); 17 > 1 
a = ^4(coshr/ — 1), t = -B(sinhr/ — 77); 17 < 1, 
where the constants A and B are 

A = A/3nGp b a 3 \R\ 2 , B = A\R\, 
R is the curvature radius and the integrals are 
r - f + T dt B 



[ v (t + t) - v (t)] , 



and 

dt 



B 
A 2 " 
B 
A 2 " 



cot 



coth 



2 



cot 



V (t + r) 



coth 



2 

7?(t + r) 



Q > 1 
Q < 1. 



The above operators can be used to build leapfrog 
or higher order symplectic integrators for particle 
motion in comoving coordinates. 



4 An energy can still be denned, as in the Layzer-Irvine 
energy equation. See Peebles, 1980. 
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